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ABSTRACT 


The numerical simulations of oscillating plus mean flow past a circular cylinder have 
‘been carried out in detail through the use of a commercially available software produced 
by CFDRC, running on a Silicon Graphics Inc. Indigo 2 Extreme computer. The 
Reynolds number, Keulegan-Carpenter number, and relative current velocity were 
systematically varied. Sensitivity analysis was performed to delineate the effects of time 
step, turbulence model and numerical schemes. The results have been compared to those 
obtained experimentally and to those predicted by the Morison Equation. In many cases 


the predicted force coefficients have shown good agreement with those obtained 


experimentally. 
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I. INTRODUCTION 


The study of time-dependent flows past bluff bodies has historically been the focus 
of'a great deal of scientific attention owing to its relevance to many and diverse applications, 
ranging from submerged structures to hot wire anemometers. Some forty years ago, 
Keulegan and Carpenter commenced the methodical investigation of the fluid-structure 
interactions which occur when bodies are immersed in unsteadily flowing fluids. Today, the 
effect of the parameters relevant to the problem, such as Reynolds number, Keulegan- 
Carpenter number and relative roughness, to name a few, is much better understood thanks 
to ongoing research in this area. In spite of this progress, real time prediction of the forces 
caused by unsteady flows on submerged objects, such as those acting on an underwater 
robotic arm or an offshore oil rig, seriously challenges our current theoretical and 
computational capabilities. 

One of the principal empirical tools used by the engineer to solve the problems 
described above, Morison's equation, is only reliable in predicting forces and moments in 
highly idealized conditions with very low or very high flow oscillation periods. 
Additionally, most real-life problems involve flows in the turbulent regime, further 
increasing the level of difficulty of both analytical and numerical solutions. 

Given the cost and complexity of the laboratory apparatus required to reproduce real 
life flow conditions, theoretical and experimental advances have been paralleled by a 
considerable research effort in the area of computational fluid dynamics to assist in the 
solution of problems related to bodies immersed in unsteadily flowing fluids. Progress in 
numerical techniques and the ever increasing power of present day computers are finally 
making it possible to overcome the barrier historically imposed by the physical and 
numerical instabilities which have caused the modelling of turbulent flows particularly 
challenging in the past. 

Whereas, until recently, numerical experimentation has only been successful in 


determining flow patterns, Strouhal numbers and force coefficients for flow regimes within 





selected ranges of relatively low Reynolds numbers, it is now possible to obtain useful data 
for more turbulent flows exhibiting extensive separation. Furthermore, commercial software 
is becoming available which affords researchers much greater flexibility than the ad hoc 
codes previously generated to solve very specific categories of problems, allowing for the 
analysis of a much broader gamut of flow conditions. 

Unfortunately, in spite of such great flexibility, modern software does not yet absolve 
the user from having to fine-tune the code by the judicious selection of parameters, 
numerical techniques and turbulence models. Even then, experimentally obtained results 
cannot be exactly duplicated in all cases. 

If complete and accurate solutions are not yet always achievable, however, 
approximate solutions can certainly aid in expanding our current level of understanding and 
provide the engineer with a valuable tool to improve design optimization. From the above 
discussion, it can be inferred that a benchmark body of reliable calculations, performed using 
the most appropriate set of computation tools, is required to calibrate our current empirical 
equations and experimental results. 

The objective of this investigation is to lay the foundations for such an endeavor and 
improve our current ability to predict the effects of time-dependent turbulent flows over 
circular cylinders through the use of a state of the art commercial code generated by CFD 
Research Corporation (CFDRC). Whenever possible, the results achieved have been 
compared to those obtained experimentally. Besides their intrinsic value, these results will 
aid in pointing out some of the strengths and weaknesses of the code and hopefully add to 
our understanding of the physics of flow types not previously analyzed. 








Il. BACKGROUND STUDIES 


About twenty years after Keulegan and Carpenter published their ground breaking 
work, a comprehensive series of experiments concerning sinusoidally oscillating flow about 
smooth and rough bodies was performed at the Naval Postgraduate School (Sarpkaya, 1976). 
This endeavor produced a wealth of practical results and introduced the parameter B (= 
Re/K) to assess the influence of scale in periodic flows. Among its findings, the dependence 
of the force-transfer coefficients on Re, K, and k/D was shown graphically and amply 
discussed. As a continuation to this research effort, cases were analyzed which involved a 
coexisting flow as specified by U=U, - U,, cos (21t/T), in which U, represents the steady 
mean velocity and U,, the amplitude of sinusoidal oscillations (Storm 1984). Since the mid- 
eighties, several numerical predictions of the Strouhal number, the pressure distribution, and 
the evolution of the lift and drag forces in steady and time-dependent ambient flows have 
been performed in an attempt to obtain results consistent with experimental data. Here, only 
the more recent and relevant investigations will be briefly reviewed. 

A finite-difference analysis of the Navier-Stokes equations for a sinusoidally 
oscillating ambient flow about a circular cylinder at K = 5 (Re = 1000) and K = 7 (Re = 700) 
has been attempted by Baba & Miyata (1987), assuming a physically unrealistic symmetric 
wake in both simulations. The results have shown that the calculations could be carried out 
only for short times (less than two cycles of flow oscillation) with a nonsuper computer. 

A similar method was used to analyze three cases (K = 5, 7, and 10) at Reynolds 
numbers around 10* (Murashige et al. 1989). The flow was perturbed by artificial means to 
trigger asymmetry. At K = 10, a transverse vortex street appeared, in agreement with flow 
visualization experiments. 

The transverse vortex street observed at K = 12 was also reproduced correctly and 
for the first time by Mostafa (1987) using multi-discrete vortices. However, the calculated 
forces were somewhat larger than those measured. 


In later years, numerical experiments with co-existing flows produced extremely 








interesting flow features. For relative current velocities V, (=U,/U,,) in the range of 0.7-0.8, 
the vortices shed nearly symmetrically at each cycle and formed a three-row vortex street in 
the range of V, = 0.6-0.7. For V, larger than about one, the vortex wake returned to the 
asymmetric mode, as in a regular vortex street. Although the scope of this work was limited 
to relatively low Reynolds numbers, the calculations of resistance in co-existing flows 
showed that both the inertia and the drag coefficients for K = 4-6 were in reasonable 
agreement with experimental data (Sarpkaya et al. 1992). 

Square cylinders have also been the subject of numerical experimentation. Earlier 
attempts were flawed by central differencing at large cell Reynolds numbers, which led to 
spatial oscillations ahead of the rectangle. Improvements were made by using time 
differencing and third-order upwinding numerical schemes on the convective terms. The 
technique just described was tried for Reynolds numbers under 3,000 (Davis & Moore 1982). 
In 1993, Kato and Launder used modified K-epsilon (k-e) and K-omega (k-Q) turbulence 
models on square cylinders, further improving the numerical results and obtaining a 
markedly superior behavior in the near-field region. 

Aerodynamic research on airfoils further demonstrated the potential of the k-e model. 
A simulation developed by Rogers (1994), although requiring more grid points than previous 
ones (such as the Baldwin-Barth model), was significantly better at computing maximum lift 
conditions and flap boundary-layer separation. 

Whereas most of the early numerical codes were generated by higher learning 
institutions or government research agencies (such as NASA), a number of very versatile 
computational fluid dynamics (CFD) software, suitable for design and analysis purposes, 
is currently being produced by the private sector. The code utilized for this investigation is 
a general-purpose CFD code with multi-domain solution capability issued by CFDRC. 

Initial results in applications germane to the subject of the current investigation 
indicated that the program was capable of predicting Strouhal numbers and forces for 
uniform flows at higher Reynolds numbers but questioned its ability to capture the high 


turbulence intensity levels present in the near-wake region (Singhal & Avva 1994). The 











program was further tested in uniform flow using the k-e and renormalization group (RNG) 
models at Reynolds numbers over 10°. Although it was able to predict Strouhal numbers 
reasonably well, the program was not tested for force coefficients under these flow 
conditions (Habchi & Hufford 1995). 

In the current investigation, the CFDRC software has been used to predict the forces 
acting on circular cylinders immersed in various types of time-dependent turbulent flows. 
It is felt that the data presented herein provide a good indication of the code's performance 


in this application. 











Til. NUMERICAL REPRESENTATION 


A. COMPUTATIONAL METHOD 


The CFD-ACE code simulates fluid flow by solving the partial differential equations 
(PDE) which govern the transport of flow quantities. Since the CFD-ACE theory manual 
discusses the solution methodology at length, only a brief outline of this topic will be 
presented here. 

The governing equations for the flow conditions investigated are the continuity and 
Navier-Stokes (NS) equations. For turbulent flows, the code applies Favre (density) 
averaging to the NS equations, so that each flow quantity @ is decomposed into a mean and 


a fluctuating part according to: 


o-+'; where: = [1] 


Note that overbar denotes Reynolds (time) averaging, while tilde denotes Favre averaging. 


Applying the Favre averaging procedure to the continuity and NS equations yields: 


dp oO 
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respectively. By applying the generalized Boussinesq eddy viscosity concept, the Reynolds 
stresses can be treated as a linear function of the mean strain rate, so the Favre averaged NS 


(FANS) equation can be expressed as: 
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in which n, represents the turbulent eddy viscosity and k is half the trace of the Reynolds 
stress tensor. 

In general, except for the continuity, each governing equation can be expressed in a 
common form comprised of a transient term, a convective term, a diffusive term and a 
source, not all of which will be present at all times. 

The equations corresponding to the flow quantities being analyzed are discretized by 
the code over finite control volumes and numerically integrated, thereby yielding a set of 
finite difference equations (FDE). Although algebraic, these FDEs are not linear and are 
therefore solved using an iterative process. 

As mentioned above, there is no governing PDE available for pressure. The task of 
solving for this flow quantity while satisfying the continuity equation is accomplished bya 
variation of the "Semi-Implicit Method for Pressure-Linked Equations, Consistent" or 
SIMPLEC algorithm. 

Briefly, the SIMPLEC procedure can be summarized as follows: 


1. Estimate the pressure field. 


2. Solve for the remaining flow quantities and check whether the continuity equation 
is satisfied. Since this will not generally be the case, an improved solution will be 
required. 


3. Find correction factors for velocities and density. These can be used to find a 
correction factor for pressure. 








4. Using the correction factor obtained in the preceding step, re-evaluate the pressure 
and velocity fields. 


5. Solve the discretized equations for other flow variables, including turbulence 
quantities, as required by the specific problem. 


6. Repeat steps 2 through 5 until the solution converges. 


The aforementioned variation used in the CFD-ACE code consists of repeating steps 
3 and 4 a number of times, thus updating the pressure, velocity and density fields prior to 
proceeding. These intermediate continuity or pressure-correction iterations have been found 


to enhance overall convergence for most flow problems. 


B. GRID GENERATION 


The grid used to analyze most models is shown in Figure 1. It was developed using 
the CFD-GEOM program and is a 1 meter by 1 meter square subdivided into four domains. 
A circle of diameter 0.05 meters, located at the geometric center of the grid, represents a 
cylinder of infinite length about which the fluid flows. The internal domain boundaries are 
located along the diagonals of the square. Sixty-five lines project radially from each of the 
four quadrants of the circle and meet the sides of the square at equispaced mesh points. 
Intersecting the radial lines, one-hundred and thirty other lines join mesh points located 
along the diagonals and spaced such that the grid is much finer in the vicinity of the circle 
than in the outer regions. 

Several requirements were considered during the grid selection process. First, the 
grid must be sufficiently extensive to prevent any vortices, still strong enough to influence 
the flow, from migrating beyond the boundaries. Also, while it was recognized that 
increasing the number of cells (as well as shortening the time interval between iterations) 


would result in a more accurate vorticity computation, efficiency considerations dictated a 








finer mesh only in the vicinity of the cylinder walls, where the effect of vorticity is 
strongest. Unfortunately, the vortices that separate and are swept away diffuse rapidly as 
the mesh coarsens, resulting in returning vortices which are probably weaker in the model 
than they would be in actuality. 

When the number of grid cells was quadrupled, marginal improvements in the results 
were obtained at the cost of a four-fold increase in CPU time. It was felt, therefore, that the 
grid just described represented the best compromise between accuracy and computation 


efficiency. 


c. BOUNDARY CONDITIONS 


In the CFD-ACE program, boundary conditions can be specified via drop-down 
menus or by editing the file in which all flow parameters are stored. Since this file has an 
"in" suffix, it will hereinafter be referred to as in-file. 

Four types of boundary conditions are used in all the simulations performed. The top 
and bottom sides of the flow region are designated as symmetry boundaries. This condition 
implies a zero normal velocity and a vanishing gradient normal to the boundary for all 
variables. 

Each cylinder quadrant is designated as a wall, implying, again, a zero normal 
velocity along the boundary. 

The vertical sides of the square, referred to as east and west boundaries, are both 
designated as exits with prescribed pressure, resulting in a paradoxical flow region with two 
outlets and no inlet. This apparent anomaly, present in all time-dependent flow simulations 
performed, can be explained by pointing out a peculiarity of the code nomenclature: while 
boundaries designated as inlets do not allow bidirectional flow, exits do. Several other inlet 
and exit conditions may be specified to accomodate supersonic, compressible and other types 
of flow. 

It should be noted that the flow conditions need not be the same along the entire 


length of each boundary. It is possible, for instance, to specify a higher flow velocity at 
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the northern half of an inlet or exit and a lower one along the southern half. This technique 
was found particularly useful when flow perturbations were desired to precipitate the 
asymmetrical shedding of vortices in uniform flow. 

A fourth type of boundary, called interface, is utilized to designate lines along which 


one domain meets with another within the flow region. 


D. FLOW CONTROL 


In the CFD-ACE program, velocity at an inlet or exit boundary can be specified 
either as a constant or as one of a limited number of time-dependent functions. Pull-down 
menus can be used to specify constant flows while time-dependent flows must be specified 
by editing the in-file. 


In the simulations performed, it was desired for the ambient velocity to vary as: 
2 
Ot) = U, - ae [5] 


in which U, represents a steady component, i.e., a current, and U,, the magnitude of the 
fluctuating component. Unfortunately, velocity was noted not to follow harmonic functions 
whenever such a function was invoked. CFDRC personnel acknowledged this anomaly 
although its underlying causes were not discussed. 

Flow can also be controlled by making use of the following relationship between 


pressure and velocity: 


9p 


27 
a ae U(t) [6] 


and by prescribing a time-dependent pressure at the exits. It can be shown that, if velocity 
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is to vary cosinusoidally, the amplitude of the differential pressure should be: 


Ap, = —EpU,Ax [7] 


where Ax is half the length of the flow region. The amplitudes specified at the east and west 
exits should be equal in magnitude and opposite in sign. Since pressure and velocity are 
offsct by /2 radians by virtue of Equation [6], differential pressure should vary sinusoidally, 
resultiag i a pressure which varies linearly in space while its slope varies harmonically in 
time. 

The technique just described allows the user to calculate the magnitude of the 
oscillating differential pressure required to generate the desired flow velocity. Whenever a 
current is present, however, a steady differential pressure is also required. It was empirically 


determined that this steady component is defined as: 


C 





Ap, - = p U,V, C ; where J, - = [3] 


in which C is dependent on K and on the time step used. The manner in which this constant 
was evaluated, as well as some recommended values for different K ranges and a time step 


of 0.02 seconds are discussed in Appendix B. 


E. DATA REDUCTION AND CALCULATION OF FORCE COEFFICIENTS 


Once the program algorithm has converged on a set of numerical values for the flow 
quantities sought at each time step, the results are recorded in a file whose suffix is ".out" 
and will therefore be referred to as the out-file. 

In order to reduce the data to a useful format, it was found necessary to transfer the 


numerical results pertinent to the investigation from the out-file to a new file. Since the 


12 





code calculates shear and pressure forces separately at each boundary designated as a wall, 
it is helpful if in the course of data transfer the contributions from each wall were 
consolidated into a single value for each cartesian axis. 

Once the total horizontal and vertical forces are available, in-line force and lift 


coefficients are obtained as follows: 


il [9] 





and: 


LF 
C,.- —— [10] 
0.5pU"D 


It is also possible for the code to record a given flow quantity at a specified point and 
tabulate the results in a monitor file. Both flow quantity and location to be monitored must 
be specified by the user. This capability was used extensively to record velocities in order 
to verify and later fine-tune the flow control technique described in an earlier section. 

Once the desired coefficient or flow quantity is tabulated, its trend can be plotted and 


further analyzed using one of several commercially available codes. 


F. INPUT OPTIONS 


The CFD-ACE Command Language Manual describes the basic techniques used to 
specify the code inputs and the other computational options available. Although drop-down 
menus can be used in most cases, certain types of inputs can only be specified by editing the 
in-file. This section will concentrate on specific aspects of input definition which were — 


found particularly relevant and in some cases critical to the type of flow conditions 
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investigated. A sample in-file with useful comment is presented in Appendix C. 

A very important aspect of flow control is the proper definition of the coefficients in 
the harmonic functions which determine the pressure differential across the flow region. Of 
these coefficients, C1 represents the steady state component (Ap, in Equation [8]), C2 the 
amplitude of the sinusoidally oscillating component ( Ap,, | in Equation [7]), and C3 its 
frequency in radians per second. The remaining coefficients, C4 and C5, were set to Zero as 
they are used for cosinusoidal components, not present in any of the cases examined. 

The user must specify an initial and a final time for all time-dependent flows. The 
time increment is indirectly specified by inputting the number of time steps to be used and 
should be chosen taking into account the fineness of the grid mesh and the flow velocity. 
In general it 1s not desireable to have the fluid travel a distance greater than half of a cell 
during a single time step. Luckily, a modicum of experimentation can quickly reveal 
whether the improvement in output quality is worth the longer computation time required 
when the ime step is shortened. | 

It should be noted that it may be desirable for the initial time to differ from zero 
because the program has the capability of restarting from a set of initial conditions recorded 
at the end of a previous run, allowing for two sets of output data to be spliced together and 
plotted continuously versus time. 

The solution control part of the in-file is one where the available options must be 
chosen judiciously and some trial and error may be required. In most cases, the default 
upwind spacial differencing scheme and Euler temporal scheme will work satisfactorily: 
Some preliminary testing of other options led to the observations which follow. 

Use of the Crank-Nicholson scheme cannot be prescribed using the appropriate drop- 
down menu as this will result in a syntax error within the in-file. In general, however, this 
technique does not produce appreciable changes with respect to the results achieved using 
the default Euler mehod. | 

While the upwind spacial differencing method is very robust, central dierencing 


may become unstable, especially when used in conjunction with fine grid meshes. This 
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problem can usually be overcome by the appropriate selection of a blending factor: in most 
models in which central differencing was used, a blending factor of 0.3 was specified; 0.5 
was successfully substituted in models which diverged. The greater the blending factor, the 
more the techinque will have an upwind character. | 

Other spacial differencing schemes did not work satisfactorily on the type of 
computer used. A later version of the program corrected this. © 

The number of solution iterations and sub-iterations for each flow quantity can only | 
be selected using some experimentation or previous experience and it will differ for each 
type of problem analyzed. In general, all iterative processes should be shortened as much 
as possible without compromising accuracy. This can only be done by solving a problem 
with a known solution and changing the number of iterations until a satisfactory result is 
obtained. The program allows for visualization of the trend in the residuals of each flow 
quantity calculation, s0 the slope of residual graphs is also helpul in determining the 
optimum number of iterations. 

| The code offers a choice of six different turbulence models, only three of which were 

utilized in the course of this investigation. The standard k-e model was found to work 
satisfactorily in most circumstances. Very minor changes were noted in the output when the 
RNG model was used. Whenever the grid-flow velocity combination is such that y’ (=(y/v 
)(t /p )°°) is smaller than about 11 at the grid cell center closest to the cylinder wall, use of 
the Low Reynolds model is appropriate and results in better numerical accuracy. Since in 
time-dependent flows velocity varies, the user needs to determine whether the values of y* 
warrants selecting the Low Reynolds model. All the above models require the user to input 
two parameters, turbulence kinetic energy and turbulence length scale. It was found that 
varying the former from a value corresponding to 3% to 10% turbulence caused minimal 
changes in the results, so 3% was used throughout. The latter was varied between 0.09 D/2 
and 0.09 D (D being the cylinder diameter) without noting significant differences. 
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IV. DISCUSSION OF RESULTS 


A. INTRODUCTION 


The results will be presented in the following order. First, a discussion of the grid 
and other code inputs used in the calculations. Second, a detailed discussion of the in-line 
and transverse force coefficients as functions of the selected governing parameters. Finally, 
a comparison, whenever possible, of the calculated and experimental values. This section 
_ is then completed with plots of vorticity and streamlines for representative values of flow 
parameters. Throughout, a special attempt is made to concentrate on the most informative 


range of the calculations rather than on a massive presentation of the data. 


B. GRID AND INPUT PARAMETERS 


As discussed in Section III B, the code used allowed one to generate various types 
of grids. For the problem under consideration, the axial symmetry of the cylinder, the 
problem's two-dimensionality and oscillating nature (left and right boundaries), and the 
extreme importance of the accurate calculation of the velocity gradients at and near the 
cylinder surface were kept in mind. Several trials and tribulations finally resulted in the grid 
shown in Fig. 1. It consists of four domains. The smallest cell size is 3.05*10° D along the 
circumference of the cylinder and 2.56*10° D in the radial direction. The largest cell, at the 
outer boundaries, is 0.313 D by 0.245 D. The size of the computational domain is 20 D by 
20 D. 

Careful consideration had to be given to the selection of the grid size for reasons 
beyond the accurate calculation of the velocity gradient on the cylinder. The fact that the 
flow oscillates makes the vortices generated during one half-cycle return toward their 
inception point. Thus, if a vortex does not retain its integrity as it transits areas of coarser 
mesh, then it will return to the cylinder with artificial diffusion above and beyond that 
imposed by the prevailing turbulence. Furthermore, the stability of the flow, regardless of 


where the cell may be, is dictated, among other parameters, by the Von Neumann criterion. 
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It states that flow at a given point should not displace more than one half of the cell on which 
the fluid particle resides. This is a difficult condition to satisfy even in unidirectional steady 
flows. Since in periodic flows, with or without current, the magnitude of the ambient 
velocity fluctuates between zero and a prescribed maximum, then one either allows violation 
of the Von Neumann criterion at times of peak velocity, so that global instability is not 
allowed to dominate the flow, or is penalized in terms of grid fineness, and consequently 
CPU time, by satisfying the criterion throughout. The experience of others, as well as the 
present one have demonstrated that the latter is impractical. 

In the present simulations, a compromise was made to maintain sufficient accuracy 
while sacrificing neither physics nor a reasonable CPU time. It is as a result of the foregoing 
considerations that the quantity U,,At/C, where C is the grid cell's characteristic length, was 
allowed to reach values on the order of 50 in the vicinity of the cylinder walls. It should be 
noted, of course, that the quantity u(t)At/C, i.e. the ratio of the actual distance traveled by the 
fluid to the computational cell's characteristic length, is expected to be substantially less than 
U,,AU/G in the immediate vicinity of the solid wall boundary. 

In the results to be presented, the flow parameters chosen were the Keulegan- 
Carpenter number K, the relative velocity parameter V, and the frequency parameter 6 or 
the Reynolds number Re. The majority of the calculations were performed with a time 
interval At of 0.02 seconds. As will be noted in connection with the sensitivity analysis 
discussion, the time interval, and other code and problem related parameters, such as the 
turbulence models, were judiciously varied. It is fully realized that a numerical simulation, 
like a physical experiment, produces only a data point on the basis of the input parameters 
chosen and discretization schemes employed in the code. Thus, the insight that one develops 
through familiarity with the phenomenon is brought to bear on the selection of the input 
parameter and on the comparison of the results with those obtained experimentally. While 
the calculations undoubtedly suffered from boundary constraints and artificial diffusion, the 
experiments have their own corresponding limitations with more or less uncertainties. 


Once the parameters were chosen, the code was capable of producing the in-line and 
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transverse forces and the necessary information for the plot of vorticity and streamlines. The 


calculations have been carried out for a minimum of seven and a maximum of 30 cycles. 


C. IN-LINE AND TRANSVERSE FORCE COEFFICIENTS 


These will be discussed for B = 1772, 2487, and 4122. In case of B = 1772, the 
analysis was carried out for six V, values, corresponding to six K values ranging from 8.6 
to 36. For B = 2487, for five V, values corresponding to five K values ranging from 10.5 to 
46. Finally, for B = 4122, for three values of V,, corresponding to three values of K ranging 
from 5.2 to 24. In addition, for B = 1772 only, the parameter V, (= KV, = U,T/D) was 
assigned three different values in order to compare the results with those obtained 
experimentally wherever possible. For the other two B values, only single values of V, were 
used. In fact, it is because of this reason that there was a one-to-one correspondence between 
V,and K for each combination of B and V,. 

Figures 2a through 2d show the in-line force coefficient (C, ) and the transverse force 
coefficient (C,) for B=1772, K=8.6 and V,=2.05. The basic difference between the first two 
and the second two figures is that the former were produced using the k-e turbulence model 
and the latter using the low Reynolds number (Low-Re) model. The purpose of this exercise 
was to demonstrate the consequences of the use of two turbulence models. 

It is realized that, for the case being discussed, the maximum Re is approximately 
15,000 and that, during any given cycle, the flow undergoes a transition from laminar to 
turbulent. Thus, the use of a turbulence model such as the k-e, tailored for use in high Re 
steady flows, may not be entirely appropriate for an unsteady flow at relatively low Re. A 
comparison of Figs. 2a and 2c shows that, whereas in the former C, does not have an 
inflection point even as late as the fourth cycle, in the latter an inflection point is evident 
shortly after the first cycle. The fact that the evolution of the inflection point is directly 
related to the occurrence of asymmetric vortex shedding is an indication of the evolution of 
the lift force. In fact, a comparison of Figs.2b and 2d shows that there is no significant lift 


in the former relative to the latter. In other words, the use of the Low-Re model, thought to 
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be more appropriate for low Re flows, gives rise to transverse force and vortex shedding after 
a few cycles from the inception of flow. 

A comparison of the turbulence models cannot be made for all the cases considered 
in this study because of obvious time limitations. Nevertheless, an effort is made to combine 
such sensitivity analyses into a section to be presented later: 

_ Returning to the presentation of the results, and in particular to the case of 8 = 1772, 
Figs. 3 and 4 show the C;, for V,=2.05 for K=15 and 36.5. It is evident from a brief perusal 
of these figures that the in-line force assumes a quasi-steady state within four to six cycles. 
The shape of the C, graph, particularly in Fig. 4, is reminiscent of an odd harmonic function, 
with the inflection points near the upper left and the lower right of each cycle associated, as 
noted earlier, with the alternate shedding of the vortices from the top and bottom of the 
cylinder. 

Figures 5a, 6a, 6c and 7 show C;, again for B=1772 for V,=4 and for three different 
K values. The most noteworthy among the four plots is Fig. 6a, where C,, exhibits irregular 
behavior after approximately eight cycles, having earlier reached a quasi-steady state. This 
seemingly strange occurrence is due to the phenomenon known as the transverse vortex 
street. For K values in the range of 10<k<13, the vortices shed from the cylinder move in 
the transverse direction and only on one side of the cylinder. This means that the vortices 
shed in each cycle go around the cylinder and find their way in the positive or negative 
vertical direction, thereby creating asymmetrical and unsteady in-line forces. This 
phenomenon has been discussed in detail by Sarpkaya and Isaacson (Sarpkaya and Isaacson, 
1981). In fact, Fig. 6b shows that a transverse force develops after approximately seven 
cycles and reflects the shedding of the vortices at the top and bottom of the cylinder for a 
number of cycles. It should be noted that Figs. 6a and 6b were developed using the 
renormalization group theory (RNG) turbulence model, whereas the low-Re model was used 
to generate Figs. 6c and 6d. Even though the C,'s do not differ significantly, the transverse 
force coefficients (C,,) show the dramatic effect of the choice of turbulence model used in the 


evolution of the vortex asymmetry. As will be noted later, not all predictions of different 
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turbulence models will yield the same result even when the same set of governing parameters 
is used. 

Figures 8 and 9a show C;,, again for B=1772, for a higher value of V, and for K=10.3 
and 36.5. Figure 9a is particularly interesting in that C, is asymmetrical in shape and 
magnitude even though K is as large as 36.5. Normally, one would expect a steady flow like 
behavior for V,=0 and relatively large K values. An examination of Figs. 9a and 9b indicates 
that neither C;, nor C, are symmetrical in the presence of a current, whereas Figs. 9c and 94d, 
which correspond to flows with no current but otherwise at identical values of K and 8, show 
that both C, and C, are as symmetrical as they can be after an initial stage of flow 
development. Thus, the presence of current has a significant impact on the topology of 
vortex shedding although the superposition of collinear waves and current do not necessarily 
lead to periodically alternating vortex shedding. It is because of this reason that C,, in Fig. 
9a exhibits symmetry neither in magnitude nor in shape. 

The next case to be considered is for an intermediate B value of 2487. Figures 10 
through 12 show C;, for V,=5.01 and for three values of K. The most significant aspect of 
the graphs is that the amplitude of the oscillation decreases with increasing K and the shape 
of the C;, trace becomes less sinusoidal, reflecting, at each half cycle, the asymmetric 
shedding of the vortices. 

Figure 13 shows the C;, for the highest B value dealt with in this investigation for 
V\=2.07. The force oscillations are, as expected, fairly sinusoidal since, for small K values, 
one does not expect significant vortex shedding. 

Figures 14a and 14b for V,=2.07 and B=4122 shows what happens to in-line and lift 
forces in the most dramatic range of K. When K is in the range between 10 and 20, vortex 
shedding becomes very chaotic, C, decreases significantly and C,, drops down relative to the 
case of K=5. 

The presence of current further complicates the flow topology with respect to the no- 
current case in this most sensitive region of the K values. The flow does not have any sort 


of symmetry about either axis. Even though the prediction of forces in this region is rather 
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difficult, there 1s one beneficial effect deriving from the smallness and the chaotic 
character of the lift, i.e. the minimization or practical elimination of the tendency of the 
cylinder to develop a dynamic response or hydroelastic oscillations. However, for K=24 
(B=4122 and V,=2.07), the vortex shedding becomes quasi-periodic and the C, increases 
considerably as shown in Figs. 15a and 15b. 

A harmonic analysis of the lift coefficient has shown that the ratio of the 
frequency of vortex shedding to the flow oscillation frequency is in the range of 3 to 4 
depending on V.,.. Such an analysis is not shown here but may be found in Storm's work. 


(Storm, 1984). 


D. FLOW VISUALIZATION 


Figures 16a through 16f show the evolution of vorticity for K=8.78 and B=1772. 
' In this particular flow regime, vortices are formed and shed from the body in quasi- 
symmetrical pairs, so that no significant amount of lift is generated. This observation is 
reinforced by examining Figs. 17a through 17f, which show the streamlines corresponding 
to the Sane flow regime. 

Figures 18a through 18f show the more complex flow topology resulting from a 
slightly higher K and a much higher B. (10.5 and 2487 respectively). The corresponding 
streamlines, shown in Figs. 19a through 19f, clearly indicate that the vortex pattern is 


asymmetrical, so that a sizable lift force results. 


E. NUMERICAL RESULTS 


The values of the maximum in-line force coefficients and the corresponding 
governing parameters are listed in Tables 1 through 5. Unfortunately, experimental results 
were available only for one set of the governing parameters. The general trends noted are 
that the agreement with both the Morison's equation and with experimental data is better 
for K values less than about 10 and larger then 20. This confirms that in a K range close 
to that in which a transverse vortex street can be expected to develop, the flow topology 


is particularly challenging to capture. Another trend noted was that, as the value of V, 
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increases, the agreement between the computed C,; values and those produced by the 
Morison's equation will deteriorate slightly. The availability of additional experimental data 


at higher V, would be required to draw more insight from the observed trend. 


F. SENSITIVITY ANALYSIS 


A fairly detailed sensitivity analysis for no-current cases was conducted by Hanson 
(Hanson, 1995). Since in that investigation the same computational grid, turbulent 
parameters, and time increment were used as those employed in this study, most of the 
findings are valid for the cases examined herein. 

The analysis conducted within the scope of this work involves primarily turbulent 
models and time increments. Since, due to time constraints, this analysis was performed on 
representative cases only, some of the observations made will be mostly qualitative in nature. 

In all cases examined, the turbulence parameters prescribed corresponded to 3% 
turbulence and an integral turbulent length scale of 0.00225. The fluid density was 998 
Kg/m? and the kinematic viscosity 1*10° m?/s. A spacial central differencing scheme with 
a blending factor of 0.3 was used, along with an Euler temporal scheme. Other 
computational parameters can be obtained from Appendix C, which is a fairly typical in-file. 

Several turbulence models were tested for the case corresponding to K=12.19, 
B=1772, V=4. While these all produced values of C,, within 3% of each other, the k-e 
model produced no inflection point after the second oscillation and the vortex pattern 
appeared to be still symmetrical after 7 cycles of flow oscillation. In contrast, the traces 
produced by the RNG and Low Re models are shown in Figs. 6a through 6d. 

As previously stated, the quantity U,,At/¢ was allowed to reach values in excess of 
50 at the cylinder walls. The case corresponding to K=24, B=4122 (Re=100,000), for which 
the initial value of U,,At/¢ was 62.3 was examined as it represents one of the highest values 
of that quantity encountered in this study. The time increment was first adjusted so that 
U,,At/¢ was approximately 12.5. This resulted in a reduction in C, from 1.22 to 1.05. (A 
similar procedure for K=46 and B=2487 resulted in a reduction in C, from 1.1 to 0.8.) A 
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further reduction in U,,At/¢ to 0.623 did not result in a sizeable change with respect to the 
case corresponding to U,,At/¢=12.5 (less than 1%). Very little effect on C,, (<2%) was also 
noted when U,,At/¢ was changed from approximately 15 to less than 0.5 for a case where 
K=12.19 and B=1772. 

It is reasonably safe, therefore, to conclude that, for the particular computational grid 
used, values of U,,At/c on the order of 15 should be acceptable. Although there is not 
sufficient data to pinpoint the critical value of U,,At/¢ for which a significant impact on the 


result should be expected, this value will almost certainly be between about 20 and 50. 
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TABLE 1 


Computed in-line force coefficients for different values of K and 6 = 1772, 
V, = UK = 2.05 


C;, (Morison) 
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TABLE 2 


Computed and experimental in-line force coefficients for different values of K and B=1772, 
V, = U,K = 4.0 





TABLE 3 


Computed in-line force coefficients for different values of K and Bp = 1772, 
V, = U,K = 6.05 
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TABLE 4 


Computed in-line force coefficients for different values of K and 8 = 2487, 
V, = U,K =5.01 





TABLE 5 


Computed in-line force coefficients for different values of K and B = 4122, 
V,. = UK = 2.07 





Zl 
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V. CONCLUSIONS 


The purpose of this investigation was a critical assessment of the numerical analyses 


of two-dimensional, sinusoidally oscillating, turbulent flows at relatively large Re, K and V; 


A number of turbulence models, including the standard k-e, RNG based k-e, and Low-Re 


models have been employed. The emphasis was on the validation and accuracy of the time- 


dependent computations. Observations relating to both numerical and physical validation 


lead to the following remarks: 


(1) 


(2) 


(3) 


(4) 


The finite differencing formulation of the Favre Averaged Navier Stokes 
equation (FANS) can reasonably solve flows within a wide range of Re using a 
modified SIMPLEC method. However, the predicted forces will almost always 
be smaller than those obtained experimentally. 


The turbulence models used do not fully capture vortex strength and prematurely 
dissipate the vortices relative to experimental observations and measurements. 
All models used herein (standard k-e, RNG based k-e, and Low-Re models) fail 
to predict accurately the exact size of the vortices and the high turbulence 
intensity levels present in the near-wake of time-dependent flows (subjected to 
unsteady pressure gradients). This is the primary cause of the aforementioned 
tendency to underpredict the force coefficients. 


The use of much finer grids, reduced time increments and higher-order spacial 
and temporal differencing schemes (at the expense of increased CPU time) will 
not always increase the accuracy of the predictions since the validity of the 
turbulence models for flows subjected to extra strains such as the time- 
dependent pressure gradients remains unknown. 


Turbulence is at present, and is likely to remain for a long time, the greatest 
roadblock to computational fluid dynamics. Some appreciation of the 
incomplete knowledge bases (both numerical and experimental), retrofitting of 
data, an the assessment of their consequences are necessary to achieve often 
qualitative and occasionally quantitative simulations, particularly in time- 
dependent flows. This is a compromise between expectations and achievables 
and between physically relevant dynamics and specific quantitative results. 
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Figure 2a. In-line force coefficient versus time. K=8.6; B=1772; V,=2.05 
(k-e turbulence model used) 
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Figure 2b. Transverse force coefficient versus time. K=8.6; B=1772; V,=2.05 
(k-€ turbulence model used) 
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Figure 2c. In-line force coefficient versus time. K=8.6; B=1772; V,=2.05 
(Low Re turbulence model used) 
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Figure 2d. Transverse force coefficient versus time. K=8.6; B=1772; V,=2.05 
(Low Re turbulence model used) 
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Figure 3. In-line force coefficient versus time. K=15; B=1772; V,=2.05 
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Figure 4. In-line force coefficient versus time. K=36.5; B=1772; V,=2.05 
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Figure 5a. In-line force coefficient versus time. K=8.78; B=1772; V,=4 
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Figure 5b. Transverse force coefficient versus time. K=8.78; B= 1772; V,=4 
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Figure 6a. In-line force coefficient versus time. K=12.19; B=1772; V,=4. 
Inception of asymmetric vortex shedding at t/T~7. (RNG turb. model) 
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Figure 6b. Transverse force coefficient versus time. K=12.19; B=1772; V,=4. 
Inception of asymmetric vortex shedding at t/T=7. (RNG turb. model) 
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Figure 6c. In-line force coefficient versus time. K=12.19; B=1772; V,=4. 
Inception of asymmetric vortex shedding at t/T~3. (Low Re turb. model) 
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Figure 6d. Transverse force coefficient versus time. K=12.19; B=1772; V|=4. 
Inception of asymmetric vortex shedding at t/T~3. (Low Re turb. model) 
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Figure 7. In-line force coefficient versus time. K=26.07; B=1772; V,=4 
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Figure 8. In-line force coefficient versus time. K=10.3; B=1772; V,=6.05 
Inception of asymmetric vortex shedding at t/T=4 
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Figure 9a. In-line force coefficient versus time. K=36.5; B=1772; V,=6.05 


Lift coeff 





0 1 2 3 4 5 6 7 
UT 


Figure 9b. Transverse force coefficient versus time. K=36.5; B=1772; V,=6.05 
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_ Figure 9c. In-line force coefficient versus time. K=36.5; B=1772; V,=0 
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Figure 9d. Transverse force coefficient versus time. K=36.5; B=1772: V_=0 
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Figure 10. In-line force coefficient versus time. K=10.5; B=2487; V,=5.01 
Inception of asymmetric vortex shedding at t/T=5 
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Figure 11. In-line force coefficient versus time. K=17.2; B=2487; V,=5.01 
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Figure 12. In-line force coefficient versus time. K=46; B=2487; V,=5.01 
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Figure. 13. In-line force coefficient versus time. K=5.1; B=4122; V,=2.07 
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Figure 14a. In-line force coefficient versus time. K=15; B=4122; V,=2.07 
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Figure 14b. Transverse force coefficient versus time. K=15; B=4122; V,=2.07 
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Figure 15a. In-line force coefficient versus time. K=24; B=4122; V,=2.07 
Inception of asymmetric vortex shedding at t/T=4 
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Figure 15b. Transverse force coefficient versus time. K=24; B=4122; V,=2.07 
Inception of asymmetric vortex shedding at t/T=4 
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Figure 17. Streamlines. K=8.78; B=1772; V,=0.43 
a. U(t)/Um = 0.99; b. U(t)/Um = 0.50; c. U(t)/Um = -0.02; d. U(t)/Um = -0.42; e. Ult)/Um = -0.58: £. U(t/Um = -0.46. 
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Figure 19. Streamlines. K=10.5; B=2487; V,=0.47 
a. U(tVUm = 0.11; b. U(t)/Um = 0.98; c. U(t)/Um = 1.46; d. U(tVUm = 1.22; e. U(tW/Um = 0.47; f. U(t)/Um = -0.29. 
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APPENDIX B. | 
RECOMMENDED VALUES FOR THE CURRENT VELOCITY- 
PRESSURE DIFFERENTIAL CORRELATION COEFFICIENT 
It was found experimentally that the coefficient C in equation [8] of section II D 
assumes different values as K and the computational time interval change. Since a 
considerable number of calculations were performed with a time interval of 0.02 seconds, 
the following table of recommended values for C, valid for time steps on the order of 
0.02 seconds, is included. It should be realized that the recommended values constitute 
only a starting point and that further fine tuning might be necessary depending on how 
accurately the value of U, is to be duplicated. One can infer from the table that, as K 
increases, C becomes closer to 0.5. In fact, when 0.5 was used in one steady flow case - 
(which can be thought as a flow with an infinite value of K), the correct value of ambient 
velocity was obtained. 
As far as the dependence of C on the time interval, as a thumbrule C will assume 
values closer to 0.5 as the time interval is shortened. Unfortunately, the nature of this 
dependence has been noted only toward the end of this investigation so that not even an 


approximate quantitative definition of it can be formulated at this time. 
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APPENDIX C. SAMPLE IN-FILE 


TITLE 'K=36.5, Beta=1772, C1=158.49, Vk=6.05, T=1.41 coarse grid, no disturbance' 


* 
* 
* 
* 
fe 2fe fe 2 afc fe 2h of ofc oft ofe afe fe afc ofe 2fe ofc oft fe ofc ofc of Model Hie 2 oft aie oe 2 oft oi ofc fe ofc ae ae I I fe Oe ae ofc 3h 3h oe He 
* 


* 


MODEL tcg115t 

FUNCTIONS 

PERIODIC WAVEPTOP C1=158.49 C2=2874.69 C3=4.453 C4=0 C5=0 
PERIODIC WAVEPBOT C1=158.49 C2=2874.69 C3=4.453 C4=0 C5=0 
PERIODIC WAVENTOP C1=-158.49 C2=-2874.69 C3=4.453 C4=0 C5=0 
PERIODIC WAVENBOT C1=-158.49 C2=-2874.69 C3=4.453 C4=0 C5=0 
END 


* 


GEOMETRY 
GRID 2D BFC 
READ GRID FROM tcg.PFG — 
* No Cell Types (Blockage or Solid) specified. 
END | 
* 
PROBLEM _ TYPE 
SOLVE FLOW TURBULENCE 
UNSTEADY TF = 0 TL = 10 STEPS = 500 
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END 
; 
PROPERTIES 
DENSITY CONSTANT 998 
VISCOSITY CONSTANT KINEMATIC 1e-06 
END 
; 
MODELS 
TURBULENCE_MODEL KE 
END 


oe 


*** Boundary Conditions *** 
** 
BOUNDARY CONDITIONS 
* boundary condition: INTERFACE 1 1 
INTERFACE 1 129 1 1 SOUTH 
* boundary condition: wall 
WALL 1 1 1 64 WEST 
U=0V=0 
* boundary condition: outritghttop 
EXIT P 129 129 32 64 EAST 
U=0V=0 P=FUNCTION WAVENTOP K = .000753 D =0 L = 0.00225 T = 
. 293 
* boundary condition: outrightbot 
EXIT P 129 129 1 31 EAST 
_ U=0V=0 P=FUNCTION_WAVENBOT K = .000753 D =0 L = 0.00225 T = 
293 | 
4 boundary condition: INTERFACE 1 2 
INTERFACE 1 129 64 64 NORTH 


52 








DOMAIN 2 
* boundary condition: INTERFACE 2 1 
INTERFACE 64 64 1 129 EAST 
* boundary condition: wall 
WALL 1 6411 SOUTH 
U=0V=0 
* boundary condition: symm 
SYMMETRY 1 64 129 129 NORTH 
* boundary condition: INTERFACE 2 2 
INTERFACE 1 1 1 129 WEST 
DOMAIN 3 
* boundary condition: INTERFACE 3_1 
INTERFACE 1 12911 SOUTH 
* boundary condition: outlefttop 
_ EXIT P 113264 WEST 
U=0V=0 P=FUNCTION _WAVEPTOP K = .000753 D = 0 L = 0.00225 T = 
293 : | 
* boundary condition: outleftbot 
EXIT P 11131 WEST 
U=0V=0 P=FUNCTION WAVEPBOT K = .000753 D =0L = 0.00225 T = 
293 | | | 
* boundary condition: wall 
WALL 129 129 1 64 EAST 
U=0V=0 
* boundary condition: INTERFACE 3 2 
INTERFACE 1 129 64 64 NORTH 
DOMAIN 4 | 
* boundary condition: INTERFACE_4_1 
INTERFACE 64 64 1 129 EAST 
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* boundary condition: symm 
SYMMETRY 1 6411 SOUTH 
* boundary condition: wall 
WALL 1 64 129 129 NORTH 
U=0V=0 
* boundary condition: INTERFACE 4 2 
INTERFACE 1 1 1 129 WEST 
* No Momentum Resistances Defined 
END 
x 
INITIAL_CONDITIONS 
* Full field initial conditions 
| U = -1.079 V =0 P=0 T = 293 K = .000753 D = 0 L = 0.00225 
* RESTART FROM tcg43.2000.AUR 
END 
* 
Ae 3 Ae Ae He fe fe fe 3h fe ofc fe fe 2h fe fe fe aft 9c 9k Solve aie fe fs 2 ye ae fe aft fe oe ofc of ofe fe ofc ofc of afc ae of ofc ac ok 
a 


* 


SOLUTION_CONTROL 
ALGORITHM SIMPLEC 
S SCHEME UPWIND RHO K D 
S SCHEME CENTRAL U V 
S BLENDING 0.3 U V 
T_SCHEME EULER 
ITERATIONS 10 
C_ITERATIONS 1 
SOLVER WHOLE_I U V PPK D 
S ITERATIONS 4 U V 
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S ITERATIONS 20 PP 
S ITERATIONS 10 KD 
INERTIAL FACTOR 0.2 U V 
INERTIAL FACTOR 0.2 K D 
RELAX 0.4 P 
RELAX 1 RHO T VIS 
MINVAL -le+20 U V 
MINVAL -le+20 P 
MINVAL 1e-06 RHO 
MINVAL le-10 T VIS 
MINVAL le-10 K D 
MAXVAL 1e+20 U V 
MAXVAL 1e+20 P RHO 
MAXVAL 5000 T 
MAXVAL let+20 VIS 
MAXVAL 1e+20 K D 

END 

‘ 

OUTPUT 
P_FORCE ON PMIN = -1.e+05 
PLOT3D ON FORMATTED 
SCALAR FILE 1 RHO P T K STRM 

DIAGNOSTICS OFF 
MONITOR 2 32 129 U 
MONITOR 4 321 U 
TIME_SAVE 250 
UNIQUE_NAME ON 

END 


> 
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